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ABSTRACT 


This report is about the effects of non-normality on the efficiency of the 
Kalman filter, particularly when the distribution of measurement errors is still 
symmetric but the tails are extended, which means that the observations are outlier- 


prone. 


“a, 
oN 


brry 
\ 
H Copy 
LONSE ora 
NN 





Accession For A 

NTIS GRA&I Ww 
DTIC TAB 
Unannournced 


Justificatio 







EI 



















By. 
Diatributions 


Avallabllity Codes 





THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research 
may not have been exercised for all cases of interest. While every effort has been 
made, within the time available, to ensure that the programs are free of 


computational and logic errors, they cannot be considered validated. Any application 


of these programs without additional verification is at the risk of the user. 
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I. INTRODUCTION 


Most contemporary automatic target trackers use a Kalman filter in their 
kernels. One of the main assumptions behind the Kalman filter algorithm is that the 
errors in the observations of the target are normally distributed. This report is about 
the effects of non-normality on the efficiency of the Kalman filter, particularly when 
the distribution of measurement errors is still symmetric but the tails are extended, 
which means that the observations are outlier-prone. 

Chapter II explains the bases of digital filters (in particular the Kalman filter) 
and the Integrated Ornstein-Uhlenbeck motion model, which are the tools for this 
study. Chapter III describes the three models that will be used in the experiments. A 
fourth mode) is presented to demonstrate the data fusion of two sources of 
information. Chapter IV describes the experiments and analyses the results and 
Chapter V summarize the conclusions. 

The "real world" situation that is inspiring our model can be described as follows: 
a station equipped with an Electronic Support Measure (ESM) equipment is trying to 
track a target by measuring its bearing and signal intensity . Those two measures have 
different inherent error, intrinsically characteristic of the equipment and of the 
physical conditions. The signal intensity provides the "range measurement” based on 
emission characteristics (such as transmission power, frequency, etc.) and propagation 
conditions. The conversion of signal intensity into range is certainly the largest source 


of error for the target state estimate. The bearing is measured directly in the 





equipment, and we expect it to be more accurate than range measurement; the exact 


precision is an equipment attribute (Figure 1). 





Figure 1 - Hypothetical situation 


In the situation described, large relative range errors and small bearing errors 


should be expected. 


Il. BACKGROUND 


This chapter was inserted in the thesis to help the reader who is not familiar 
with the Kalman filter and the Integrated Ornstein-Uhlenbeck model. For the reader 
who is already familiar with these subjects, this section is optional. The following is a 


synopsis of various documents [Ref. 1 through 4). 


A. DIGITAL FILTERS' 

A sequence of bearings and ranges is produced when periodic measures of a 
target are made by a Radar, Sonar or other detection device. Those observations carry 
in addition to their inaccuracy an associated uncertainty which is usually represented 
by additive noise [Ref. 5]. 

The algorithms used to process those measurements are called Digital Filters. 
If we call the inputs IN(t) and the outputs OUT(t), Equation 1 represents a digital 
filter 


OUT) = XY C, ING-) + YD, OUTL-A (2.1) 
i=0 


i=l 


Where the C,’s and D,’s are constants. The formula says that the filter is a linear 
combination of measures and previous values obtained in the filter. If all coefficients 
D, are null, the filter is called non-recursive. In this thesis we are interested in a 
recursive filter, the Kalman Filter, which uses the last measurement and a 


combination of all the previous ones. 


1 This section is an adapted translation of Reference 1. 





The name "filter" comes from the analogy with the idea of a filter in electronics, 
where for example, when we want to eliminate high frequency signals we use a "Low 
Pass" filter, and only slow changes in the signal are allowed to go through (Figure 2). 


Low pass characteristics are obtained on the filter in Equation 2.1 when D,> >C,. 


ONE 


Figure 2 - Electronic filter 


Analogously, our filter will compare the measure with the previous state of the system, 
allowing the last measure to have more or less influence on the system state change, 
according to the gain established. 
1. Recursive Filters 
As mentioned previously, a recursive filter uses the last measurement and 
a combination of the previous measures that is in the system state. As an example of 


a recursive algorithm, look at the equation 


OUT®) = (1 - ? IN(t) + ; OUTt-1) (32) 


Considering the input at t=0 as IN(0) and OUT(t)=0 for t<0, then tt. : sequence of 


values will be generated at each unitary time step: 








OUT(O) = IN(O) 


OUT) = (1 - » INQ) + + IN(O) 


i (2.3) 


1 
OUT(2) = (1 ry IN(2) z OUT(1) 
wie ey er eee 
(1 - 3) INQ) + (& - =) IND) + = INO 


Note that the value IN(0) is never forgotten and can have more or less influence in the 
present state depending on the value of K. Note also that the value of OUT(t-1) 
economically condenses all the previous information up to time t-1. 
a. Filtersa Banda py 
Consider a stationary target (order zero system) for which measures 
of one coordinate are being taken. One filter that might be used to represent the 


system is 


OUT) 


a IN) + (l-a) OUTU-L 
OUTKt-1) + a(IN(t) - OUT(t-1)) 


(2.4) 


In this case the actual state (target position) is obtained by combining the previous 
state with the difference between the measurement and the previous System state, 
multiplied by the gain a. If the gain a = 1 total confidence is placed in the measure and 
the previous state is forgotten. Conversely a =0 will not use the measure, keeping the 
previous state. This is called the « filter. 

Consider now a target with constant velocity (first order system)(Figure 


3). The filter now has to predict the position at the time ¢ of the measure before 


comparing it with the estimate. For that system, calling X the target position and V, 








the target velocity, using the same idea as for the previous filter the equations are: 


PREDICTION X,() = X,¢-An) + V, At (2.5) 
ESTIMATION XO) = XO + (XO - X,0) (2.6) 
where X, = predicted position 


X,, = measured position 
X, = estimated position 
a = filter gain 
That means the estimated position is placed somewhere between the predicted 


position and the measured position, using a linear combination. 





Figure 3 - Prediction Step. 


It is possible now, using the previous estimated position and the measured position, 
to calculate the measured velocity. Then a linear combination of the measured and 
predicted velocities is formed using a new weighing parameter, the gain B, to achieve 
a better estimate of velocity. Remember that the predicted velocity at time t is equal 


to the estimated velocity at time t-At by assumption. The equation for the estimated 


velocity is: 


VO = VO + BVO - VO) (2.7) 

or, in a continuous time formulation 

£X0 = 2x (0-a9) + Pxo ~X0) . (2.8) 
Note that when the measure is greater than prediction, the correction in the 
estimated velocity is positive (as it should be). The parameter 6 is the filter velocity 
gain and the filter is called «_f. In the same way a filter for a system with constant 
acceleration (second order system) might be derived. A new gain factor will appear and 
the filter will be called « § y [Ref. 1). 

Looking at the filter a_f, it is intuitive that the first measures will 
have a high value in the target state determination, since they are the only 
information the filter has. It is also intuitive that the choice for the gain value will 
affect the speed of convergence of the filter to the correct estimates of position and 
velocity. We now seek a way to calculate the "best" (in some sense) gain to be used in 
the filter at each measurement. That is the work performed by the Kalman filter. 

b. The Kalman Filter 

Consider a system of order zero (stationary target) from which two 
measures X,, where 1 is the time, are obtained. Knowing that the true position is X,,if 
the measures are corrupted by a Gaussian noise and the errors are independent, the 


first two measures are: 


X, =X +e,, 
(2.9) 
and X, =X, +e, 


Estimating the target position by a linear combination of the measures, the estimate 


and its error are: 


X, = aX, + bX, 
(2.10) 
e=X,-X, , 
where a and 0 are the linear coefficients. The criteria chosen to select the "best" 


estimate are minimum average error and minimum square error. In that case the 


equation for a null average error, using the facts that the mean for the errors is null 


and E[X,] = Xp) is: 


Efe] = Efax, + bX, - X)) 
Efa(X, + €,) + W(X, + €,) - Xl 


= aX, + bX, - X, = 90 


(2.11) 


a+b=] 
The equation for the average square error using the result of Equation 2.11 is: 


Efe] = El[(a+b-1)X, + ae, + be,}’] 
(2.12) 


= Ef(ae, + be,)"] as (a+b) = 1 





Which can be also written as: 


Ele’ 


Ejae? + 2ab ee, + b7e3] 


a*E[e7] + b?E[e?] (2.13) 
= a’g? + bo, = a207 + (1-a)*or . 


After differentiating and equating to zero to find the minimum, the gains and the 


minimum square error will be: 











i 
2 2 
ae OFS oj 
GO) + 0, me + a 
oO, 02 
(2.14) 
1 
o of 
b= lang ad 
a + a; ae + A 
o; 02 
2 
9192 1 
Efe’] = oo 
0; + 03 2 + i ; (2.15) 
GO; 02 


Note that the minimum square error is smaller than the variance for the two 


measures, which means the variance for the estimate is better than the variance of the 
two measures. Note also that if 
0,>06, , thn a>b, (2.16) 


which means that more weight is given to the more precise measure. These results 


will be applied for the minimum mean square error in a recursive filter that also 


combines two values, the measure and the prediction. The best gain for each step is 


given by Equation 2.14. 


Preeti aes (2.17) 


variance of the measure 


= 
g 
a 
J 
“ 


variance of the prediction 


Q 
i] 


For a first order system (constant speed) the process of estimation has 
to be preceded by a prediction (Figure 3 and Equations 5 and 6). The intent now is to 
calculate « and B following a minimum mean square error criteria of optimization that 
is 

Flex] is minimum .. @ is the best 


Eley, 7] is minimum: 8 is the best 


where €y = X, - X__ error in estimated position 


e, = Vx - Vy error in estimated velocity 


aS 
' 


true velocity 


X = true position 


At the instant of the measurement a prediction operation is performed. The values for 


the position and the error in predicted position are: 








X (t+ At) XO + Vy At 


(2.18) 
Ey = ey, + €y, At . 
The equations for the velocity and velocity error are: 
Vy (+A) = Vy. @ 
(2.19) 
&y, = &y 
The equations for the mean square errors are: 
Fle] = Ele,’ + dey ey, + €y, I 
(2.20) 
‘ Oy” + 2oyy,? + oy? 
Eley 7] = Fley,"] = oy,” : (2.21) 


But these are the values for variance of position after prediction and variance of 


velocity after prediction. 


2 2 2 
Oy, = Oy2 + 20 o 
Xp x XVx, Vy, ees; 
2 2 
o 
Vip vy, 


At time ¢ a measure is obtained and a new estimate is constructed. The equations for 


position and variance of position can be derived as [Ref. 1): 








X, = X, + a(X, - X,) 





2 202 (2.23) 
Pa Go = G Ox, 
Oy. = (1-a)oy, = (1 - —*—)o,, = —**- 
Oy, +0, Oy, + Om 


The equations for velocity and variance of velocity after estimation can be derived as 


[Ref. 1]: 
Vy 


‘e 


Fey, 


V, + BX, - X,) 
(2.24) 


2 2 242 2,2 
Crm * 2Boyy, + B*o,, + Boy, : 


After differentiating the last equation and equating it to zero, the equation for the 


velocity gain is: 


2 


oO 
fs 2 (2.25) 
0, + Ox 
So, after the operation of estimation, the variance will be 
Oy,” = Fey, 7] = Op - Boy,” (2.26) 


To complete the algorithm it is necessary to determine the covariances between 
position and velocity for the operations of prediction and estimation. Those can be 
derived after some computation as [Ref. 1]: 
Oxy, = Eley €y] = Ef(ey. 7 €y, ev, J = Ory, * Oy,” 
(2.27) 
Ory, = Flex €y,] =(1 - O)0xy, 
The equations described are used sequentially in predictions (or 


movements) and estimations (or measurements) to update the state vector and 
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covariance matrix. The system condition is then revised recursively at each time step 
either by a movement and a measurement or just by a movement when no observation 
is available. 

"All of Kalman’s computations can be thought of as manipulations of 
(multivariate) Normal probability distributions" [Ref. 6). The attractive feature of the 
linear plant model is that normality is preserved. 

The Kalman filter, in the way it was described in this section, works 
for targets with constant course and speed. If the target maneuvers after the filter has 
settled with an optimal gain, the measures will not affect the gain and the filter will 
not follow the maneuver. To prevent that “disconnection from reality” the realistic 
deviation of the target from a precisely straight, or other deterministic track is 
modeled by a random process. 

It is important not to confuse the model noise, also called plant noise, 
with the measure noise. The measure noise appears due to the inaccuracy and 
uncertainty in the measurement and is simulated by a random number generator. The 
plant noise is placed in the model to reflect reality and is achieved by the plant noise 
covariance matrix Q. 

To simplify the calculation, the Kalman filter can be represented in 
matrix notation. The operations performed are the same, but the manipulation is 
much easier. A primer for the subject is Reference 6, from where the notation will be 
borrowed. 

The process can be divided in two steps, the measurement step that 
corresponds to the prediction, and the movement step that corresponds to the 


estimation. The equations for the steps, omitting the time index, are: 
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movement step 
B= op + by (2.28) 


B= 247+ @ 


measurement step 


K = SH™(H=H?™ + Ry! 

(2.29) 
B =p + K(Z-p,-Hp) 
Z=(U-KHZU , 


where 


® is the plant matrix 
(p,.Q) is the mean and covariance of the plant noise 
H is the measurement matrix 
(p,R) is the mean and covariance of the measurement noise 
(u,L) is the mean and covariance of the state of the system 
is the measurement vector 
K is the Kalman gain matrix 


I is the identity matrix . 
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B. MOTION MODELS 

When tracking a target, detection equipment will make observations at intervals 
of time, which may or may not be equal. Between those measures, the tracker will 
have to draw inferences about the target’s position. For that, the tracker must assume 
that the target is moving in a certain way, which is called the motion model (Figure 


4). 





Figure 4 - Motion Model Examples. 


The motion model reflects the way the target is expected to behave. It forecasts 
the position, compares with the observation, and reveals the likely or unlikely track. 
The simplest motion model is the constant course and speed model, when the tracker 
assumes the target is moving in a straight line with constant velocity. But many other 
possibilities exist. For example the GST (Generic Statistical Tracker) uses four motion 
models [Ref. 4]. The GST keeps track of both the target position (in two dimensions) 
and velocity. By changing some parameters in the tracker one can obtain the motion 
models described in Figure 4. As can be deduced, target’s position prediction is as good 


as the coincidence of the motion model used and the target real motion. 
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Many motion models are not physically realizable by any ship. One model that 
is physically realizable is the random tour where the target moves at constant speed 
and executes random (uniform from 0 to 2x) course changes at exponentially 
distributed random times. The random tour can be approximated by a Gauss-Markov 
process called Ornstein-Uhlenbeck process. This approximation is actually very good, 
under certain assumptions and the process can be used as a motion model with many 
advantages [Ref. 3]. 

1. The Ornstein-Uhlenbeck Process” 

The Ornstein-Uuilenbeck is a special case of the Generalized Langevin 
Model which is a model based in linear stochastic differential equations and can be 


described by the equations: 


dX(t) = Vy(Odt 
dV (t) = -BV,(t)dt + K(x) dt + (dWit) 


(2.30) 


where X(t) = X coordinate at time t 
V(t) = X coordinate of velocity at time t 
K(x,t) = an acceleration produced by an external force field 
dWit) = white noise (Wiener differential) 
€ = scale parameter 
B = parameter . 

The Langevin model is used to describe the motion of a particle subjected 
to a frictional damping force with coefficient B>0, an external force field and random 
shocks described by a white noise. The Ornstein-Uhlenbeck process is defined when 
there is no external field (K(x,t)=0 ) and it roughly corresponds to an unconstrained 


random walk without drift. 


2 This section is an adapted transcription of Reference 2. 
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The O.U. model becomes: 


dX(t) = V(t) dt 
(2.31) 
dY(t) = Vd) dt 
dV f(t) = -B,V(pdt + C,dw() 
(2.32) 


dV {t) = -BVfOdt + C,dW,() 


where B,,B,,¢, and (, are positive constants 
dW,(t) and dW,(t) are white noise 


The O.U. process becomes I.0.U. (Integrated Ornstein-Uhlenbeck) process 
when Equations 2.31 and 2.32 are integrated to compute the values for position and 
velocity. Looking at the velocity expression (Equation 2.30) one can see that the 
velocity component satisfies a stochastic differential equation that reflects 
deterministic and random influences. The first term decelerates the target at a rate 
proportional to its velocity, B being the constant of proportionality, and is the 
deterministic part. Additionally, the target undergoes a rapidly varying random 
acceleration which is idealized as a white noise. After a time, the first term will decay 
to zero and the random term will dominate. That means the ability to predict the 
velocity at a future time will decay asymptotically to zero, with B as the parameter 
that controls the rate at which the deterministic part becomes dominated by the 
random component. Although the differential equation characterization of the O.U. 
process is not hard to solve it will not be derived in this work. The computation 
developed in this report will use a discretized version [Ref 4]. The equations for the 


1.0.U. process in matrix notation will appear later in Chapter III. 
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It is important to understand the exact meaning of the parameters (¢ and 
B of the I.0.U. motion model . The I.0.U. process involves white noise and, in the 
technical terminology of filtering theory, is non-realizable. No ship of any kind, 
merchant or warship, that has ever been built can move as required by the O.U. or 
1.0.U. stochastic differential equations models. By contrast, the Random tour family 
of models are all physically realizable by actual ships. The approximation of a R.T. by 
an 1.0.U. motion model is actually very good and a correspondence between the 
parameters can be established. The Random tour is generated by a target moving with 
constant velocity V and changing course at intervals of time exponentially distributed 
with mean time between course alterations 1/4. The courses are chosen independently 
from a uniform distribution [0,27]. Thus the average number of course changes in (0,t] 
is At. It can be shown that the correspondence between the R.T. and I.0.U. parameters 


is 


(2.33) 


so the parameters B and ¢ can be described as: 


B = effective average rate of target course 
change per unit time 


Bvt RMS target speed 


The I.0.U. process will be used in this work to provide a motion model for targets with 


random behavior as will be described. 





Ill. MODEL DESCRIPTION 


When we talk about target tracking, normally a single filter will not respond 
perfectly to all situations. We will now describe four filters for three different motion 
models. These filter: will be developed to accomplish our test purpose which is to 
examine the effects of non-normal] outlier-prone observational errors on Kalman 
filtered tracks. They are not optimal for all situations. The test will be performed with 
three tracks: Inbound constant course and speed target; Crossing constant course and 
speed target; and Random Tour target. The fourth filter will be developed to 


demonstrate the use of two sources of information tracking a single target. 


A. DESCRIPTION OF THE ONE SOURCE (SINGLE TARGET) MODELS 
1. The Range, Bearing, Radia] and Angular Velocity Model (Filter 

Model 1) 

The problem geometry is defined by the observations on range and bearing 
of a target from a stationary source, positioned at the origin of a Cartesian coordinate 
system (Figure 5). It seems natural to keep the same orientation for a model. The 
state vector for that case would have four components; range, bearing, radial and 
angular velocity and can be represented as: 


B=(Upboty by) (3.1) 
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This approach tries to keep the measurement as a linear function of the state vector 


and thus to avoid the use of an extended Kalman filter. 


Figure 5 - Problem Geometry. 


The equations for the Kalman filter in matrix notation, were described in 
Chapter II (Equations 2.28 and 2.29). Those equations will be employed to derive the 
matrix to be used in the model. This model will be used on experiment 1 (inbound 
CCS target). 
a. The movement step 
The state vector and the covariance matrix are updated at each 
movement step. First the Kalman filter "moves" the target in accordance with the 
plant model, estimating the position just before the next measurement (just like a 
dead reckoning would do using the last position course and speed). Second it computes 
the expansion that the time elapsed between steps causes in the target’s region of 


uncertainty. To find the plant model matrix (®) we look at the equations of that model: 


p(t+AD= p,() + Hy (At + noise 
peltt+Ad= p(t) + by At + noise 

(3.2) 
Hy t+AD= py (© + noise 


Hy, (t+A= py (+ noise 


These equations represent a target movement from time (t) to time 
(t+ A), using the velocity computed in the previous state and adding a plant noise that 
will be null in this model. The plant matrix has to satisfy the equations described, 
which can be summarized as p(t+ A) = ® p(t) + noise. The result for the plant matrix 
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1 0 At O 

sic (3.3) 
001 0 
0001 


The plant noise covariance matrix Q, is defined for this model as a null matrix. 
b. Measurement Step 
The measurement step will first compute the Kalman gain matrix K. 
The Kalman gain determines how seriously the new measured value will be taken and 
how extensively the estimate will "bounce"; a small K makes the filter lethargic. The 
state vector and the state covariance matrix are then updated. The first gives the best 
estimate of the target at the moment and the second shrinks the area of uncertainty. 


The state covariance matrix reflects the uncertainty of the state vector components 
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and with the state vector, condenses all the information in the system state. The 


measurement matrix H relates the 2 measurements to the state vector by the equation 


Z=HX+V (3.4) 


where V is the measurement noise. In our model H will be a 2X4 matrix of the form 


a-|? 0 0 | (3.5) 
0 1 0 0 
2. The (X, Y, V,, V,) Simple Model (Filter Model 2) 

A model using the Cartesian coordinate system has some general 
advantages over the previous model such as: the transition to a model with two sensors 
is sometimes less complicated than for one in polar coordinates; a target with constant 
course and speed has constant velocities in X and Y; the process is independent for the 
coordinates X and Y, so the computation can be easily divided for parallel processing. 
The only problem is to keep the observational variables (range and bearing) a linear 
combination of the state vector (X,Y,V, and V,). A conversion of coordinates from polar 
coordinates to Cartesian coordinates is used prior to the application of the data into 
the filter, making the filter believe that everything is in coordinates X and Y [Ref. 7]. 


The equations that relate range and bearing with the coordinates X and Y are well 


known to be: 








X = R sin(8) 
Y = R cos) . 


(3.6) 


"On account of the non-linear transformation (Equation 3.9), the measurement errors 
on the Cartesian coordinates are non-Gaussian distributed and the optimal filter would 
be non-linear" [Ref. 7]. In order to avoid this difficulty, assume as a reasonable 
hypothesis that the errors in polar coordinates are small compared with the true target 
coordinates. Under this assumption the Cartesian errors are obtained by 
differentiation of Equation 3.9 and the relationship between the errors in both 
coordinates system is linear and the Gaussian probability distribution is maintained. 
It is also necessary to convert the variances in range and bearing to variances in X and 


Y. The formulas are: 
o% = sin%(@) 0% + R? cos*(®) 06 
0% = cos%(@) 0% + R? sin%(8) 07 (3.7) 
Oxy = (on-R?06) sin® cos® . 

A model using those transformations can be easily derived in a way similar to that 


done for the first model. 


a. Movement Step 
The state vector is composed of the coordinates X and Y and the 


respective velocities. The notation used is: 


B= (hy > By ’ By, ry by)” d (3.8) 


To define the plant matrix it is necessary to look at the equations of the plant model 
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pyt+Ad) = py) + py(OAt + noise 
pyt+AD = pO + uy (DAF + noise 
(3.9) 
By(t+Ad) = py, + noise 
By (t+Ad) = Py + noise . 
A plant matrix ® equal to that of the previous model (Equation 3.3) will satisfy these 
equations. The plant noise for the purpose of the test (CCS target) will be null. 

b. Measurement Step 

The measurement covariance matrix R uses the converted variances 
from range and bearing to X and Y (Equation 3.7). As the matrix H simply extracts the 
first two components of the state vector and the conversion from the measurement in 
range and bearing to X and Y is prior to the use of the data, the matrix H is the same 
for all models (Equation 3.4). 

3. The (X, Y, V,, V,) Model Using I.0.U. Process (Filter Model 3). 

As mentioned in Chapter II, the 1.0.U. motion model may be a good 
approximation to the Random Tour motion model, still allowing one to remain within 
the class of Gaussian approximations to target motion. The modifications from the 
previous model are only in the movement step, so we will not explain the 
measurement step. 

a. Movement Step 

To derive the plant model matrix (®) for the model we first look at the 


equations that describe the motion model: 
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1 
X(t+At) = X(t) + ——(1-e B44) V, At 
(t+Ad) ®*+ 5 0-¢ \ Ve 


1 
Y(t+Ad) = YX) + ——(1-e 844 V, At 
(t+ An) © + orf ey) Vy 


(3.10) 
V,ft+Ad = eB Vin 
Vit+Ad = eB Vn 
So the plant matrix is: 
1 0O Obi 0 
>= oe eh (3.11) 
0 0 5b O 
0 0 O- b2 
where 
b1 = 4(1-e*4 
b2 = e Bar 


The plant noise to be injected in the system has mean zero and covariance matrix Q 


defined as [Ref. 4]: 


Q- 0 cd oO «2 (3.15) 
0 c2. «(0 c3 


where 
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- 1/6, ~ diy — 4,-Bat . -2B4 
el 3) [2Ar Ge 4e eo)] 
€2 = a1) (1 - eBay? (3.16) 
c3 = 5 £) (i-e B45 


A Fortran code for this model is presented in Appendix A. 


B. DESCRIPTION OF THE 2 SOURCES (SINGLE TARGET) MODEL. 
Considering the same type of sensor, two sources of information will normally 
give more information than just one. When more than one source is available and they 
are separated in space, the necessary link and information fusion may degrade the 
system reliability . Sometimes the second source is not available at all and the problem 
has to be solved for just one platform, optimizing the information available. "The 
combination of 2 sources of information can be realized in a track selection fusion, in 
which the tracker with the best solution is selected and no real fusion is computed, 
or in a state vector fusion, in which the associated state vectors are combined in a 
linear estimator to derive a central-level state estimate” (Ref. 8]. The estimator in this 
case must consider sensor-specific parameters and can be implemented in a fashion 
similar to the Kalman filter to form a minimum mean square estimate. Another 
possibility is measurement fusion that uses the measures to compute a central estimate 
of target state. In this case we have to consider if the measures are synchronous or 
not, as that changes the way we combine them before the application into the filter. 
In this work a state vector fusion using a Kalman filter with the I.0.U. motion 


model is implemented. In this case only one fused state vector and fused covariance 
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matrix is employed to represent the system. The problem geometry is shown in Figure 
6. The sources are assumed to have the same orientation (positive X axis is to the 
right). Note that the error ellipses for each individual source’s observations don’t have 
the same orientation, so the measurement covariance matrix R is different for each 
source. However, conversion of the variances from range and bearing to X and Y 
(Equation 3.10) before using the data in the filter, allows the matrix R, although 
reflecting measures with different orientations, can be used by both sources with no 
distinction to update a unique fused state vector. The good part is that the Kalman 
filter will not just combine measures with different orientation but will also take into 
account the precision reflected in the values of the matrix R and will give more weight 


to the better measure. 


source 1 measure 
true 
position 


target path 


source 2 
measure 





Figure 6 - Two sources problem geometry 
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IV - MODEL EVALUATION 


A. TARGET MOTION GENERATOR 

Evaluation of a tracker requires a target motion generator to create possible 
paths that a target could produce. Noise must be added to simulate the observed track. 
The observations arrive as a sequence of bearings and ranges of the target, so after 
computing the X and Y coordinates of each motion step, we calculate the range and 
bearing from the tracker unit and then add the noise. This is an attempt to represent 
actual measured physical conditions, specifically by assuming a range error with mean 
zero and standard deviation equal to a percentage of the current target real distance, 
and a normal bearing error with mean zero and constant standard deviation. A Fortran 
code of the motion generator is presented in Appendix B, and a more detailed 
description follows. 

1. Target Measurement Error 

The target motion generator has 2 different measurement noises that can 

be used, a Normal noise and a non-Normal, outlier-prone noise simulated by 
sculpturing a Normal distribution. The normal noise is generated using the method 
of Box-Muller [Ref. 9] . The bearing standard deviation is 0.01 radian, and that for the 


range is a constant percentage of the actual target range (input for the program). 


Those figures can, of course, be modified to any desired value. 








Asculptured noise [Ref. 10] is a modification that stretches the tails of the 
normal noise to any desired value of kurtosis*. This device models the situation when 
the noise is not normal but still symmetric. It represents the occurrence of gross 
errors or outliers. The process of distributional sculpturing, described in Reference 6, 


suggests two shaping functions for stretching symmetric distributions: 


(@) S(Z=1+hZ* h>o 
(4.1) 

(i) S(=e ho . 
If Z is a random variable normally distributed with mean 0 and standard deviation o, 
then Y = ZxS(Z) resembles the distribution of Z for small values, but lengthens the 
tails of the distribution for large Z. 

To obtain the desired kurtosis in the sculptured distribution it is necessary 

to analyze the moments (particularly the second and fourth) of that distribution. They 
are represented (Equations 4.2 and 4.3) in terms of those for the original Z for each 


transformation. 


@) Y= Z(1+hZ?) 


m,(¥) = m,(Z)+h m,(Z) = 0 
(4.2) 
m,(Y) = m,(Z)+2h m,(Z)+h?m(Z) 


mY) = m(Z)+4h m(Z)+6h?m,(Z)+4h? m,(Z)+h* m,,(Z) 


3 It is also possible to sculpture for skewness, but this is 
not done in this work. 








Gi) Y= Ze (4.3) 


mY) = _ oe he (4.4) 
: (1-4ho?)? 40? 
30° 1 
Yy = ———— h<— Fs (4.5) 
= (1-8h0?)? 80? 


Where o” is the variance of Z [Ref. 10]. Those formulas are applicable to any Normal 
distribution. If h exceeds the indicated limits (Equations 4.4 and 4.5) the moments 
become infinite; even so the models may be of interest. The moments for the Gaussian 
distribution with mean 0 and standard deviation o can be found by expanding both 


sides of the equality 


Jeo? 
me} =e? (4.6) 


and comparing the terms with same degree in 6. It can be seen that the odd moments 


are null and the even moments are given by the formula 


Ez" ~ n! o 
n o 
G 2 





(4.7) 


The tail extension can be measured by the kurtosis, obtained by dividing 
the fourth moment by the square of the second moment ,and then subtracting the 


constant factor 3. Calculation of the fourth moment for transformation (i) needs all the 


even moments up to the twelfth of the original distribution (Normal in our case). 











E[Z7] = 0” 


E[Z‘] = 304 
Z = 150° 
i (4.8) 
E{Z*] = 1050° 
E[Z"] = 9450"° 


E[Z"] = 103950” 

The purpose of this work is to compare the behavior of the tracker when 
using an observational error noise that is Normal with mean 0 and standard deviation 
o, with the performance when the error noise is sculptured with the same mean and 
standard deviation but with a given kurtosis. The procedure is to start with a normal 
distribution that, after the transformation, will be sculptured with standard deviation 
of o and a given kurtosis. An algorithm that carries this on is as follows: 

1. Generate a random variable Z that is N(0,1). 

2. Pick an appropriate value of sculpturing parameter h. 

3. Apply the transformation to the random variable Z to get Y. Divide by 
the transformed standard deviation and then multiply by the desired standard 


deviation op. 


Y 
ie gi. 4. 
21 Spy °* a 


After those steps the variance of Z, is op and the kurtosis is 


Fizz) _ 


(Z,) = 3. (4.10) 





E\Z7? 








The attractive feature of the outlined solution is that the original distribution is 
N(0,1), simplifying the calculation of the moments for the sculptured distribution. 
Applying the solution to the first transformation: 

1. Define a N(0,1) random variable Z 


2. Define a random variable W and a random variable Zr as 


w = —ZU+hz?) 
VAR{Z(1 +hZ?)] (4.11) 
Z,= Wo, 


At this point the variance of Z, is o,” and the kurtosis is: 


on AW) | EW) 








y12) = ———— -3 
"(oh EW)? ] 
(4.12) 
4 
_¢(_za+2 )) _, 
VVAR[Z(1 +hZ7)] 
After some computation the Kurtosis can be expressed in terms of h as: 
3 

11Z,) = 10395h* + 4(945)h> + 6(105)h? + 4(15)h + 3 -~3 (4.13) 


(15h? + 6h + 1) 
The appropriate solution for Equation 4.13 will give the correct value for 
h to be used in Equation 4.11. The results of solving the equation numerically for 


different values of kurtosis are presented in Table 1. 
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TABLE 1 - SOLUTIONS FOR EQUATION 4.12 


Kurtosis 





Similarly, solutions for the second transformation can be found. The values for h, for 
both transformations, were confirmed by a simulation using APL2/32. 
2. Target Motion Pattern. 

The target motion generator has 4 different motion patterns, each with a 
particular characteristic to test. All 4 patterns have the same output to the file 
MOTIOx RES (the ’x’ is the motion model number) containing 9 variables described 
below: 

- TIME - Simulation time (hours). 

- X - Target x true coordinate (nautical miles). 


- Y - Target y true coordinate (nautical miles). 





- RANGE - Target true range (nautical miles). 

- BEARI - Target true bearing (degrees with 000 being the positive x axis). 

- XN - Target x coordinate with normal error (nautical miles). 

- YN - Target y coordinate with normal error (nautical miles). 

- RANGEN - Target range with normal error (nautical miles). 

- BEARIN - Target bearing with norma) error (as in BEARI). 

Only the 2 patterns that will be used in the experiments will be described. 

a. Motion Pattern 1 - Constant Course and Speed Target 

This pattern accepts as input the target initial position (coordinates X 

and Y in nautical miles), course and speed. It then updates the target position at each 
time step by using the Equation 4.14. The result is a straight path with the course 


computed counterclockwise from the positive X axis and constant speed. 


X(t+Ad) = X(t) + Vy At 
(4.14) 
Y(t+At) = Y() + Vy At 


Where V, and Vy are the velocities in X and Y respectively. 
b. Motion Pattern 4 - Random Tour Target 
This motion model was described in Chapter II. The inputs are the 


target initial coordinates, speed, and the parameter A for the exponential distribution 


(the input is 1/4). 





B. EVALUATION. 


To gain an insight into the problem, which is to examine the effect of non- 
normal outlier-prone observational error distributions in the Kalman filter track, we 
will perform 3 experiments using a 2 level factorial design (Ref. 7]. In a 2 level factorial 
design a set of orthogonal dummy variables is created and a linear regression is 
performed. Each dummy variable has 2 possible values (+1 or -1) that correspond to 
the high and low level of the factors. The model used for the linear regression of the 


dummy variables on the response variable being analyzed is 


RESP = 9+ ,4, +74, +134, +0,4,4, +4, d, +1644, +14,4,4, 
WHERE 
RESP is the response variable, (4.15) 


d, are the dummy variables for K, SD and AT and 


N, are the regression coefficients. 


The computed regression coefficients will possibly indicate how the response variable 
will behave when we change the factors from low to high level (for further information 
see Box, Hunter and Hunter (Ref. 11]). Three factors will be used in the design: 

1. Kurtosis of the range error distribution. 

2. Standard deviation of the range error distribution. 

3. Time interval between observations. 

One experiment consists of 2* combinations of factors and levels. A certain 
number of runs (500 or 1000) is performed for each combination. Each run consists of 
a sequence of observations. For each observation, the squared distance between the 


true position and the estimated position, D’, is computed. The distances for one 
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complete run are then averaged, and that will be used as the experiment response 
variable. So, the response variable is the expected value for D? (E[D’)) for one run that 
will be called ED’. The sample variance for the response variable is also computed in 
the simulations. The levels for the factors are presented in Table 2. 


TABLE 2 - FACTORS AND LEVELS 


KURTOSIS (K) 0 (Normal) 2 (3 for the first 
experiment) 













The standard deviation for bearing errors is constant and has the value 0.01 
Radians. 
1. Experiment 1: Inbound Target 
An inbound target is running at constant course 225° (000° is the positive 
X axis) and constant speed 10 and is tracked by filter model 1. Arun starts at 60 miles 
from the source and stops at 40 miles. A number of these runs is performed until we 
reach a total of 200,000 observations for each combination. That means 500 runs of 400 
observations each for AT =.005 and 1000 runs of 200 observations each for AT=.01. The 
state vector is initialized for each run with the first observation and null velocities. 


The covariance matrix is also initialized at each run with the values: 
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The results are presented in Table 3. Note that the response variable has small values, 


considering that the mean distance of the target from the source is 50 miles. As we 
are using a null plant noise covariance matrix, after a short while the filter tracks with 
precision and gives a very small weight to the observations. 


TABLE 3 - RESULTS FOR EXPERIMENT 1 


7 ae | 
Ba ol 
Ba a 





| 
| 
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The pooled sample variance is defined as (Ref. 7]: 


2 2 2 2 
52 = TS; + T,S; + TS, + ......... +T,S; (4.17) 
Dot Fy Ty > snes +T, 
with T=T1,+T, +... +r; af (4.18) 


where T', = Gf: for the estimate of variance for set of runs i 


= # of runs - 1 


After computation the pooled variance is 0.18045. Since each main effect 
is a subtraction of two averages, with each average containing eight observations, the 


variance of each effect is: 


> + t We 
Vege = UX. - XJ = (E+ OS = 


4 (4.16) 
= 0.045 

The main effects and interactions were computed and the results are 
presented in Table 4. The values found in Table 4 can be related to the regression 
coefficients (Equation 4.15) after the inclusion of a scale factor [Ref. 11]. 

The conclusions to be extracted from this experiment are: there are no 
significant (more than 2 STD. ERROR) interaction effects, so the main effects can be 
considered separately; for the levels and variations between levels used, the standard 
deviation of the measured errors affects the track significantly (as we should expect), 
and the effect of kurtosis and interval between observations are not significant. As will 
be seen in the third experiment (1.0.U. model) the levels used affect the result of the 


experiment. 


38 





TABLE 4 - EFFECTS AND INTERACTIONS (EXPERIMENT 1) 





CT ESTIMATE STD. ‘ERROR 
pe 


2. Experiment 2: Crossing Target 
A crossing target is running at constant course 315° and constant speed 10 
and is tracked by filter model 2 ( X,Y,Vx,Vy simple model). A run starts at coordinates 
X = 26 miles and Y = 42 miles, stopping 2 hours later. A total of 500 runs of 400 
observations each for AT=0.005 and 1000 runs of 200 observations each for AT=0.01 
is performed. The state vector is initialized, for each run, with the first observation 
converted by Equation 3.6 and null velocities. The covariance matrix is also initialized 


for each run with the values: 


39 








The results are presented on Table 5. 
TABLE 5 - RESULTS FOR EXPERIMENT 2 





BRCERECE 


The response variable values are also small, but this comes as no surprise 
since this filter is also optimal for a CCS target. After computation the pooled variance 
for the response variable is 0.34 and the variance of each effect is 0.086. The main 
effects and interactions the results are presented in Table 6. The conclusions to be 


extracted from this trial are the same as for the previous experiment. A few number 
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of runs high value of ED* were observed on the eighth combination, signifying that 
the track was bad on those runs. 


TABLE 6 - EFFECTS AND INTERACTIONS (EXPERIMENT 2) 


| ESTIMATE STD. ERROR 





We repeated the experiment and counted the number of runs with ED’ 
greater than 7.5 (heuristically assumed as a threshold to consider that the track was 
lost) and none of the runs in all sets of runs reached that value. 

3. Experiment 3: Random Tour Target 
a. Single Source Model 
A target performing a Random Tour starts movement at coordinates 
X=39 and Y=39 with constant speed 10 and is tracked by filter model 3 ( XY,V,,V, 
I.0.U. model). A set of 1000 runs with 200 observations each is performed for each 


combination of factors and levels. The state vector and covariance matrix are initialized 
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as in experiment 2. The number of runs with ED? greater than 7.5 are computed. The 
results showed that for combination number 3 (line 3 of Table 7) the tracker lost track 
4 times (out of 1000). The track was also lost for combination 4 (42 times), for 
combination 7 (7 times) and for combination 8 (80 times). The variances for 
combinations 4 and 8 (consequently the pooled variance) were so big that significant 
results could not be extracted. We eliminated those runs where the track was lost and 
performed the experiment again until 1000 good runs were reached‘. The results are 
presented in Table 7. 

TABLE 7 - RESULTS FOR EXPERIMENT 3 (ONE SOURCE) 


= 





* Although this procedure is not easily statistically 
interpretable, it was used for the sake of comparison with the 
other experiments (bad runs wiil be considered later). 
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The pooled variance for the response variable is 0.53 and the variance 
of each effect is 0.1315. Results of computation of the main effects and interactions 
are presented in Table 8. The standard deviation of range error is the only effect that 
is significant, at the levels and level variations used. 


TABLE 8 - EFFECTS AND INTERACTIONS (EXPERIMENT 3, ONE SOURCE) 





The experiments gave insight into the problem. We discovered that the 
tracker is losing track in a number of runs for the sets of runs that have high levels 
of K and SD. To quantify the percentage of runs in which the track is lost, we decided 
to run two more simulations; one at 5% and other at 10% of SD. The program will now 
vary the kurtosis from 0 to 5, counting the number of tracks lost in 3 categories. The 
interval between observations is fixed at 0.01 and the bad runs will be considered in 


the computations. For 5% SD no single run lost track for all values of kurtosis. 
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For 10% SD the results are summarized in the Table 9. The fields LOST 1, LOST 2 


and LOST 38 are the number of runs ,out of 1000 runs, with ED* between 7.5 and 10, 
10 and 50 and greater than 50, respectively. The field EDN? is the average squared 


difference between the measured and the true positions. 


TABLE 9 - EFFECT OF KURTOSIS VARIATION (ONE SOURCE) 










3137 | 3141 | 3145 | 3148 | 31.50 | 3153 
VARIED*] saa beso coe 1059 | 1196 | 501.58 | 1088.95 ae 


As we increase kurtosis the following happens: 

1. The values for ED? are not increasing much (for small values of SD 
they even decrease), and in the worst case we are still improving when considering the 
estimated position compared with the measured position (in average values). 

2. What is not apparent in the average values presented is that we will 
lose track in some runs, and that the number of tracks lost increases rapidly with the 


kurtosis. 
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3. The variances are growing very fast as the number of tracks lost 
increase. This is the reason why we could not get a reasonable pooled variance to 
compare with the mean effects when the levels of kurtosis and SD were high. 

4. The values for EDN? are increasing very slowly. 

b. Two Sources Model 

We now want to access the impact of non-Normal errors on the model 
that uses two sources of information by observing how the results obtained in the 
previous section will change. A simulation similar to the last one performed on the 
previous section was executed using the 2 sources model described in Chapter III. 
Another source was positioned at coordinates X = 78 and Y = 0, and alternate 
observations are made (each source performs an observation at .01 hour interval). 
Observations from both sources had the same range error distribution. The results are 
presented in Table 10. The inclusion of an additional source of information does not 
prevent the tracker from loosing track, but there is a large reduction on the number 


of tracks lost and a considerable improvement on the track quality. 
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TABLE 10 - EFFECTS OF KURTOSIS VARIATION (TWO SOURCES) 











V - CONCLUSIONS AND RECOMMENDATIONS. 


The conclusions presented in this section reflect the results of the experiments 
performed in Chapter IV. The values obtained and the conclusions extracted are for 
the models and levels used. 

There is no significant change in the response variable when we switch from the 
lower to the higher level of kurtosis, keeping constant the other factors. For small 
values of range error standard deviation (less than 5%), the effect of tail extension in 
the error distribution is not noticeable in all three models. For larger values of range 
error standard deviation, as we increase the value of kurtosis the trackers will start to 
lose track (threshold of 7.5) at the percentages presented on Table 11 for filter model 
3. In an attempt to visualize the reasons why the tracks were lost, the seeds for the 
bad runs were recorded and some of the tracks were reproduced and plotted in 
GRAFSTAT. A number of plots showed that the first observation was an outlier and 
the tracker took a long while to recognize this fact. No other apparent reason for track 


loss (other than bad observations caused by the long tail error distribution) was found. 
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TABLE 11 - PERCENTAGE OF TRACKS LOST 


KURTOSIS PERCENTAGE OF 
RUNS LOST 





A. SUMMARY 

Kalman filtering is an extremely broad subject. The models presented can be 
embellished to account for multiple targets, false echoes and clutter, targets with 
random acceleration and many other issues that appear in the real world . The two 
(and multiple) sources model is another area for research. Comparison between 
different fusion methodologies: track selection fusion, state vector fusion, and 
measurement fusion models to find the best solution for an issue is certainly a 


challerging and practical area for further study. 





APPENDIX A - FORTRAN PROGRAM FOR FILTER MODEL 3 


PROGRAM KALXY3 
EREREEEERREREESESEEESESEDEEEEEREEEREEEEEEEESEREREEEREREEERE HSER EREEEESS 
* THIS PROGRAM SIMULATES THE OPERATION OF A KALMAN FILTER 
* USING THE COORDINATES X AND Y AND THE I.0.U. PROCESS 
* AS A MOTION MODEL. 


SHEKEKSESESEKETETAESSSSKKSECEKSESE KESSEL SKRESCSSSERESESAKSEKKSESKESKHSEKEKERESEREE 
* 


PEELE EEE EE ES DESCRIPTION FOR SOME VARIABLES SERSEEESEKEKESSEKEEEEEERE 
*BETA PARAMETER IN THE IOU PROCESS 

* B1,B2,C VARIABLES DEFINED IN THE IOU PROCESS 

* C1,C2,C3 VARIABLES DEFINED IN THE IOU PROCESS 

*DELTAT INTERVAL TIME BETWEEN OBSERVATIONS (HOURS) 

*I,J | AUXILIARY VARIABLE 

*LL,CC AUXILIARY VARIABLE 

*ID IDENTITY MATRIX (4X4) 

*K KALMAN GAIN MATRIX (4X2) 

*MU STATE ESTIMATE (4X1) 

*R COVARIANCE MATRIX OF THE MEASUREMENT NOISE (2X2) 
*SIGMA PARAMETER IN THE IOU PROCESS 

*SIG COVARIANCE MATRIX ESTIMATE (4X4) 

*SIZE NUMBER OF POINTS TO BE GENERATED IN THE SIMULATION 
*TIME SIMULATION TIME (HOURS) 

*Z OBSERVATION (2X1) 

* 


SEEKERS SRSEKSKESKKSESSSESRASKSSASSEKSSKSSSESKSSSEKKEREKHESEKKESKKESSESKSKERSEKKAESE 
* 


INTEGER SIZE,I,J,LL,CC 
REAL BETA,B1,B2,C1,C2,C3,C 
REAL DDAR(2,2),DELTAT,H(2,4),HT(4,2), HSIGHT(2,2), HMU(2,1) 
REAL AUXSIG(4,4),Q(4,4),ID(4,4), IDMIKH(4,4),SCX,SCY 
REAL SIGMA,SIGHT(4,2),SIG(4,4),SIGPLUS(4,4),SCOREX,SCOREY 
REAL K(4,2),KH(4,4),KZMIHMU(4,1), RANGE, BEARI 
REAL MU(4,1),MUPLUS(4,1),PERC,PHI(4,4),PHIT(4,4), PHISIG(4,4) 
REAL TIME2,XN, YN,Z(2,1),ZMIHMU(2,1),R(2,2),RADD(2,2) 
* INITIALIZATION 
PARAMETER (PI = 3.141593) 
DATA H/1.,0.,0.,1.,4°0./, ID/1.,4*0.,1.,4%0.,1.,4*0.,1./ 
SIZE = 200 
DELTAT = .01 
SIGMA = 10. * 2. **.5 
BETA = 2. 
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PERC = .30 
DO 120 LL=1,4 
DO 130 CC=1,4 
PHI(LL,CC) = 0. 
QULL,CC) = 0. 
IF (LL .EQ. CC) THEN 
PHI(LL,CC) = 1. 
ENDIF 
130 CONTINUE 
120 CONTINUE 
PHI(3,3) = EXP(-BETA * DELTAT) 
PHI(4,4) = PHIG,3) 
PHI(1,3) = (1 - PHI(3,3)) / BETA 
PHI(2,4) = PHI(1,3) 
Q(1,1) = .5 * (SIGMA / BETA)**2) * (2.*DELTAT - (3. - 4. * 
Cc PHI(3,3) + EXP(-2*BETA*DELTAT)) / BETA)) 
Q(2,2) = QU) 
Q(1,3) = .5 * (SIGMA/BETA) * (1-PHI(3,3))) **2 
Q(2,4) = Q(1,3) 
Q(3,D) = Q(1,3) 
Q(4,2) = Q(1,3) 
Q(3,3) = .5 * (SIGMA**2 / BETA) * (1- EXP(-2* BETA * DELTAT)) 
Q(4,4) = Q(3,3) 
OPEN(UNIT = 12,FILE="D:\APL2\MOTIO4.RES’,STATUS =’OLD’) 
* OPEN(UNIT=11,FILE=’D:\APL2\KALXY3.RES’,STATUS =’NEW’) 
MUG,1) = 0. 
MU(4,)) = 0. 
SCOREX = 0. 
SCX = 0. 
SCOREY = 0. 
SCY = 0. 
DO 170 LL=1,4 
DO 180 CC=1,4 
SIGCLL,CC) = 0. 
180 CONTINUE 
170 CONTINUE 
SIG(1,1) = 100. 
SIG(2,2) = 100. 
SIG(3,3) = 50. 
SIG(4,4) = 50. 
DO 200 I=1,SIZE 
SHSTSSESSEARSESESESSSEKESESSESSE MEASUREMENT STEP SHEESEKAKSEKSSEAHSSTKSSKSEKSERE SEE 
* — WRITE(1,") 
*  WRITE(11,*)"MEASUREMENT STEP ‘I 
READ(12,110) TIME2,X,Y,XN, YN,RANGE,BEARI 
110 FORMAT(1X,F5.2,2(1X,F6.2),19X,2(1X,F8.4),1X,F10.6,1X,F 11.6) 
BEARI = BEARI * PI /180. 
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*** TRANSFORMATION OF COORDINATES 
Z(1,1) = RANGE * COS(BEARI) 
Z(2,1) = RANGE * SIN(BEARI) 
R(1,1) = (SIN(BEARI)* PERC * RANGE)**2 + 
Cc (RANGE * COS(BEARI) * .01)**2 
R(2,2) = (COS(BEARI)* PERC * RANGE)**2 + 


Cc (RANGE * SIN(BEARD) * .01)**2 
R(1,2) = (PERC* RANGE )**2 - (RANGE* .01)**2) 
Cc * SIN(BEARI) * COS(BEARI) 


R(2,1) = R(1,2) 
*** INITIAL STATE VECTOR 
IF(LEQ.1) THEN 
MU(1,1) = Z(1,1) 
MU(2,1) = Z(2,1) 
ENDIF 
*** COMPUTATION OF THE KALMAN GAIN "K" 
CALL TRANSP(H,2,4,HT) 
CALL MATMUX(SIGHT,SIG,4,4,HT,4,2) 
CALL MATMUX(HSIGHT,H.2,4,SIGHT,4,2) 
CALL MATADD(RADD,HSIGHT,R,2,2) 
CALL INV2X2(RADD,DDAR) 
CALL MATMUX(K,SIGHT,4,2,DDAR,2,2) 
* WRITE(11,*)’K’ 
. WRITE(11,*)((K(LL,CC),CC = 1,2), LL = 1,4) 
*** COMPUTATION OF THE STATE VECTOR "MUPLUS" 
CALL MATMUX(HMU,H,2,4,MU,4,1) 
CALL MATSUB(ZMIHMU,Z,HMU,2,1) 
CALL MATMUX(KZMIHMU,K,4,2,ZMIHMU,2,1) 
CALL MATADD(MUPLUS,MU,KZMIHMU,4,1) 
: WRITE(11,*) MUPLUS’ 
- WRITE(11,*)((MUPLUS(LL,CC),CC = 1,1),LL= 1,4) 
*** COMPUTATION OF THE COVARIANCE MATRIX "SIGPLUS" 
CALL MATMUX(KH,K,4,2,H,2,4) 
CALL MATSUB(IDMIKH,ID,KH,4,4) 
CALL MATMUX(SIGPLUS,IDMIKH,4,4,SIG,4,4) 
‘ WRITE(11,*)SIGPLUS’ 
i WRITE(11,*)((SIGPLUS(LL,CC),CC = 1,4),LL = 1,4) 
SRHSHHRESRHESREKS MOVEMENT STEP FSESKSSSKEKAAKSESESEKESEEESKEEKSESEKES KCREKREEE 
. WRITE(11,*) 
. WRITE(11,*) MOVEMENT STEP’ 
*** COMPUTATION OF THE STATE VECTOR "MU" 
CALL MATMUX(MU,PHI,4,4,MUPLUS,4,1) 
* WRITE(11,*)"MU’ 
: WRITE(11,*)((MU(LL,CC),CC = 1,1),LL = 1,4) 
*** COMPUTATION OF THE COVARIANCE MATRIX "SIG" 
CALL TRANSP(PHI,4,4,PHIT) 
CALL MATMUX(PHISIG,PHI,4,4,SIGPLUS,4,4) 
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CALL MATMUX(AUXSIG,PHISIG,4,4,PHIT,4,4) 
CALL MATADD(SIG,AUXSIG,Q,4,4) 
WRITE(11,*)’SIG? 
WRITE(11,*)((SIG(LL,CC),CC= 1,4), LL=1,4) 
XE = MU(1,1) 
YE = MU(2,1) 
WRITE(11,890)X,Y,XN,YN,XE, YE 
FORMAT(6(1X,F10.4)) 
SCOREX = SCOREX + (X-XE)**2 /SIZE 
SCX = SCX + (K-XN)**2 /SIZE 
SCOREY = SCOREY + (Y-YE)**2 /SIZE 
SCY = SCY + (Y-YN)**2 /SIZE 
200 CONTINUE 

PRINT*,SCOREX,SCX,SCOREY,SCY 

STOP 

END 


* 
* 
SUBROUTINE MATMUX(RESUL,MULTP,LINP,COLP,MULTR,LINR,COLR) 
* SUBROUTINE FOR MATRIX MULTIPLICATION 
INTEGER LINP,COLP,LINR,COLR,K,L,M 
REAL MULTP(LINP,COLP),MULTR(LINR,COLR),RESUL(LINP,COLR),SOMA 
DO 120 K= 1,LINP 
DO 130 L= 1,COLR 
SOMA = 0. 
DO 140 M=1,COLP 
SOMA = SOMA + MULTP(K,M) * MULTR(M,L) 
140 CONTINUE 
RESUL(K,L) = SOMA 
130 CONTINUE 
120 CONTINUE 
RETURN 
END 


SUBROUTINE MATADD(RES,ADD1,ADD2,NROW,NCOL) 
*SUBROUTINE FOR MATRIX ADDITION 
INTEGER NROW,NCOL,K,L 
REAL RES(NROW,NCOL),ADD1(NROW,NCOL),ADD2(NROW,NCOL) 
DO 150 K=1,NROW 
DO 160 L=1,NCOL 
RES(K,L) = ADD1(K,L) + ADD2(K,L) 
160 CONTINUE 
150 CONTINUE 
RETURN 
END 


* 








SUBROUTINE MATSUB(RES,SUB1,SUB2,NROW,NCOL) 
* SUBROUTINE FOR MATRIX SUBTRACTION 
INTEGER NROW,NCOL,K,L 
REAL RES(NROW,NCOL),SUB1(NROW,NCOL),SUB2(NROW,NCOL) 
DO 170 K=1,NROW 
DO 180 L=1,NCOL 
RES(K,L) = SUB1(K,L) - SUB2(K,L) 
180 CONTINUE 
170 CONTINUE 
RETURN 
END 
* 
SUBROUTINE TRANSP(A,AROW,ACOL,T) 
* SUBROUTINE FOR MATRIX TRANSPOSITION 
INTEGER AROW,ACOL,K,L 
REAL A(AROW,ACOL),T(ACOL,AROW) 
DO 190 K=1,AROW 
DO 195 L=1,ACOL 
T(L,K) = A(K,L) 
195 CONTINUE 
190 CONTINUE 
RETURN 
END 


* 
* 


SUBROUTINE INV2X2(MAT,MATINV) 
* SUBROUTINE FOR 2X2 MATRIX INVERTION 
REAL MAT(2,2), MATINV(2,2),DET 
DET = (MAT(1,1) * MAT(2,2)) - (MAT(1,2) * MAT(2,)) 
IF(DET .EQ. 0.) THEN 
PRINT*,’MATRIX IS NOT INVERTIBLE’ 
STOP 
ENDIF 
MATINV(1,1) = MAT(2,2) / DET 
MATINV(1,2) = - MAT(1,2) / DET 
MATINV(2,1) = - MAT(2,1) / DET 
MATINV(2,2) = MAT(1,1) / DET 
RETURN 
END 
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APPENDIX B - FORTRAN PROGRAM FOR THE MOTION GENERATOR 


PROGRAM MOTION 
REESASESHESKSKESSSSSSEKEKSSESEKSSERSESSEESKESSESKESSSSSEKKESESKSERSEKEKEKESEKESREEEESE 
* THIS PROGRAM SIMULATES 4 TARGET MOTION PATTERNS , USING THE 
* NORMAL RANDOM GENERATOR OF THE PCSIMUTIL PACKAGE FROM PROF. 
* MIKE BAILEY. THE RANDOM NUMBER GENERATOR IS FOR A 32 BITS PC. 
* THE OUTPUT FOR FILE E:\MOTION.RES , WHERE THE ’N’ IS 
* SUBSTITUTED BY THE TARGET PATTERN, SHOWS 9 VARIABLES IN THE 
* FOLLOWING ORDER : 
* 1)TIME 2)X 3)Y 4)RANGE 5)BEARI 6)XN 7)YN 8)RANGEN 9)BEARIN 
* 


SHERESEEERESEESEEESESS ESSER SEKESHSEERATSERESERERESEKESEKHAESKEEKEEKRERKERERRERES 
* 
FHEESSSSEHKSERSESAEHEKEHSESESES VARIABLES SREKEKSESEKESSESEKESEHKESEKASERKKEEKERREREEEE 


*ALTT RANDOM TIME TO NEXT COURSE ALTERATION ON MODEL 4 
* ALTTIME TIME OF NEXT COURSE ALTERATION ON MODEL 4 (HOURS) 
* ALPHA 1/MEAN TIME BETWEEN COURSE ALTERATION ON MODEL 4 
* BCOURSE TARGET BASE COURSE (DEGREES) 

* BSPEED TARGET BASE SPEED (KNOTS) 

* BEARI BEARING FROM OBSERVER TO TARGET (RADIANS) 

* BEARIN BEARING FROM OBSERVER TO TARGET WITH ERROR NOISE 
if (RADIANS) 

* CHOICE TARGET MOTION MODEL CHOICE (1 TO 4) 

* COURSE TARGET COURSE (RADIANS) 

*CTIME CORRECTED TIME TO BE USED AS AN ANGLE IN RADIANS 

* DELTAT INTERVAL TIME BETWEEN OBSERVATIONS (HOURS) 
*DSEED SEED FOR RANDOM NUMBER GENERATOR 

*H SCULPTURING PARAMETER 

“1 AUXILIARY VARIABLE 

*MBCD MAXIMUM BASE COURSE DEVIATION 

*MBSD MAXIMUM BASE SPEED DEVIATION 

* MEASDEV VECTOR OF MEASURES STANDARD DEVIATION (NM,RAD) 

* NOICE NOISE ERROR CHOICE (1 OR 2) 

*NORM VECTOR WITH NORMAL RANDOM VARIABLES 

*PERC PERCENTAGE OF ACTUAL RANGE USED TO COMPUTE RANGE S.D. 
*PI ASITSELF 

* RANGE TARGET RANGE (NAUTICAL MILES) 

* RANGEN TARGET RANGE WITH ERROR NOISE (NAUTICAL MILES) 
*SPESD STANDARD DEVIATION FOR SPEED NOISE (KNOTS) 

*SPEED TARGET SPEED (KNOTS) 

*SPEEDX TARGET SPEED IN THE X AXIS (KNOTS) 

* SPEEDY TARGET SPEED IN THE Y AXIS (KNOTS) 
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*SIZE NUMBER OF POINTS TO BE GENERATED IN THE SIMULATION 
*TIME SIMULATION TIME (HOURS) 

*U AUXILIARY VARIABLE 

*X TARGET X TRUE COORDINATE (NAUTICAL MILES) 

*XN TARGET X COORDINATE WITH ERROR NOISE (NAUTICAL MILES) 
7y TARGET Y TRUE COORDINATE (NAUTICAL MILES) 

*YN TARGET Y COORDINATE WITH ERROR NOISE (NAUTICAL MILES) 
* 


SSSCSEEESSESESEESEREESSESSESSSESETSESESESESSSSESSAE SEEKS SSESEREEKSESSEREAHESEE 


* INITIALIZATION 
INTEGER SIZE,I,CHOICE,NOICE 
DOUBLE PRECISION DSEED 
REAL ALTT,ALTTIME,ALPHA 
REAL BCOURSE,BSPEED,BEARI,BEARIN,COURSE,CTIME,DELTAT,H 
REAL MBCD,MBSD,MEASDV(2),NORM(2),PI,PERC,SPEED,SPEEDX 
REAL SPEEDY,SPESD,RANGE, TIME,X.Y,XN, YN,U(1) 
PARAMETER (PI = 3.141593) 
DATA DELTAT/.01/ , SIZE /200/, H /0.1097466/ 
TIME = 0. 
ALTTIME = 0. _ 
DSEED = 2332. 
* USER CHOICE OF TARGET PATTERN 
WRITE(6,*)’ TARGET MOTION GENERATOR?’ 
WRITE@,*) 
10 WRITE(6,*)’PLEASE INPUT THE DESIRED MEASUREMENT NOISE * 
WRITE(6,*). 1- NORMAL NOISE’ 
WRITE(6,*) 2- SCULPTURED NOISE’ 


READ(6,*) NOICE 
IF((NOICE.LT.1).OR.(NOICE.GT.2)) GOTO 10 
WRITE(66,*) 
20 WRITE(6,*)’PLEASE INPUT THE DESIRED MOTION PATTERN NUMBER ’ 
WRITE(6,*)’ 1 CONSTANT COURSE AND SPEED’ 
WRITE(6,*)’ 2 SINUSOIDAL TRANSITOR TARGET ’ 
WRITE(£6,*)’ 3 SECOND ORDER GAUSS MARKOV TARGET ’ 
WRITE(66,*)’ 4 RANDOwl TOUR TARGET’ 
READ(6,*) CHOICE 


IF((CHOICE.LT.1).OR.(CHOICE.GT.4)) GOTO 20 
GOTO (100,200,300,400),CHOICE 


SESSEEESEES GENERATION OF A MODEL 1 TARGET SHERKSESSEREREASESHEHERHESS EEE EE 
® 


100 WRITE(6,*Y INPUT TARGET INITIAL X-Y POSITION (E.G., 0. 0.)’ 
READ(6,*) X,Y 
WRITE(6,*) INPUT TARGET BASE COURSE AND SPEED ’ 
READ(6,*) COURSE, SPEED 
COURSE = COURSE * PI / 180. 
SPEEDX = SPEED * COS(COURSE) 
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SPEEDY = SPEED * SIN(COURSE) 
WRITE(6,*) INPUT RANGE STANDARD DEVIATION, PERCENTAGE OF 
WRITE(6,*)) | ACTUAL RANGE (E.G. 5% IS 5) 
READ(6,*) PERC 
PERC = PERC / 100. 
OPEN(UNIT =08, FILE =’D:\APL2\TEST.RES’ , STATUS =’NEW’) 
* MAIN LOOP FOR MODEL 1 
DO 190 I=1,SIZE 
X = X + SPEEDX * DELTAT 
Y = Y + SPEEDY * DELTAT 
TIME = TIME + DELTAT 
RANGE = SQRT(X**2 + Y**2) 
BEARI = ATAN2(¥,X) 
* GENERATION OF 2 NORMAL RANDOM BEARING AND RANGE 
CALL LNORPC(DSEED,NORM,2) 
MEASDV(1) = PERC * RANGE 
MEASDV(2) = .01 
IF(NOICE.EQ.1) THEN 
RANGEN = RANGE + NORM(1) * MEASDV(1) 
BEARIN = BEARI + NORM(2) * MEASDV(2) 
ELSE 
RANGEN = RANGE + (NORM(1) * (1. + H * NGRM(1) ** 2)) 
/ (1. + 6*H + 15* H**2)**.5)) * MEASDV(1) 
BEARIN = BEARI + (NORM(2) * (1. + H * NORM(2) ** 2)) 
Cc / (GQ. + 6*H + 15* H**2)**.5)) * MEASDV(2) 
ENDIF 
XN = RANGEN * COS(BEARIN) 
YN = RANGEN * SIN(BEARIN) 
WRITE(08,115) TIME,X,Y,RANGE, 180.*BEARI/PI.XN,YN,RANGEN, 
Cc 180.*BEARIN/PI 
115 FORMAT(1X,F5.2,3(1X,F6.2),1X,F11.6,2(1X,F8.4), 
Cc 1X,F10.6,1X,F11.6) 
190 CONTINUE 
GOTO 900 


* 


SKEKKKERVEKK ERE GENERATION OF A MODEL 2 TARGET SREEKKKESEEKREKREREHEEKEREEE 
* 
200 WRITE(6,*) INPUT TARGET INITIAL X-Y POSITION (E.G., 0. 0.) 
READ(6,*) X,Y 
WRITE(6,*) INPUT TARGET BASE COURSE AND SPEED ’ 
READ(6,*) BCOURSE,BSPEED 
BCOURSE = BCOURSE * PI / 180. 
WRITE(€6,*)’ INPUT MAXIMUM DEVIATION FROM BASE COURSE (E.G., + /- 
C 45 IS 45) 
READ(6,*) MBCD 
MBCD = MBCD * PI / 180. 
WRITE(6,*) INPUT MAXIMUM DEVIATION FROM BASE SPEED (E.G., + /- 2 
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C KTSIS2/ 
READ(6,*) MBSD 
WRITE(6,*) INPUT RANGE STANDARD DEVIATION, PERCENTAGE OF ’ 
WRITE(6,*)” ACTUAL RANGE (©.G. 5% IS 5) 
READ(6,*) PERC 
PERC = PERC / 100. 
OPEN(UNIT=09, FILE =’D:\APL2\MOTIO2.RES’ , STATUS=’NEW’) 
DO 290 I = 1, SIZE 
CTIME = AMOD(TIME,2*PI) 
COURSE = BCOURSE + MBCD * SIN(SPEED*CTIME) 
SPEED = BSPEED + MBSD * SIN(CTIME) 
X = X + SPEED * COS(COURSE) * DELTAT 
Y = Y + SPEED * SIN(COURSE) * DELTAT 
RANGE = SQRT(X**2 + Y**2) 
BEARI = ATAN2(Y,X) 
* GENERATION OF 2 NORMAL RANDOM BEARING AND RANGE 
CALL LNORPC(DSEED,NORM,2) 
MEASDV(1) = PERC * RANGE 
MEASDV(2) = . 
IF(NOICE.EQ.1) THEN 
RANGEN = RANGE + NORM(1) * MEASDV(1) 
BEARIN = BEARI + NORM(2) * MEASDV(2) 
ELSE 
RANGEN = RANGE + ((NORM(1) * (1. + H * NORM(1) ** 2)) 
C / (1. + 6*H + 15* H**2)**.5)) * MEASDV(1) 
BEARIN = BEARI + ((NORM(2) * (1. + H * NORM(2) ** 2)) 
/ (1. + 6*H + 15* H**2)**.5)) * MEASDV(2) 
ENDIF 
XN = RANGEN * COS(BEARIN) 
YN = RANGEN * SIN(BEARIN) 
WRITE(09,115) TIME,X,Y,RANGE, 180.*BEARI/PILXN,YN,RANGEN, 
Cc 180.*BEARIN/PI 
TIME = TIME + DELTAT 
290 CONTINUE 
GOTO 900 


SRSSKSSEHERKRESES GENERATION OF A MODEL 3 TARGET SRERSSSEEEKESEKEREEKEREEE 
* 
300 WRITE(6,*) INPUT TARGET INITIAL X-Y POSITION (£.G., 0. 0.)’ 

READ(6,*) X,Y 
WRITE(6,*) INPUT TARGET INITIAL COURSE (DEGREES) AND SPEED (KTS)’ 
READ(6,*) BCOURSE,BSPEED 
BCOURSE = BCOURSE * PI / 180. 
WRITE(€6,*) INPUT SPEED STANDARD DEVIATION FOR NORMAL NOISE ’ 
READ(6,*) SPESD 
WRITE(6,*) INPUT RANGE STANDARD DEVIATION, PERCENTAGE OF ’ 
WRITE(6,*)’ ACTUAL RANGE (E.G. 5% IS 5) 
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READ(6,*) PERC 
PERC = PERC / 100. 
SPEEDX = BSPEED * COS(BCOURSE) 
SPEEDY = BSPEED * SIN(BCOURSE) 
OPEN(UNIT= 10, FILE =’D:\APL2\MOTIO3.RES’ , STATUS=’NEW’) 
DO 390 I = 1,SIZE 
CALL LNORPC(DSEED,NORM,z2) 
SPEEDX = .8 * SPEEDX + NORM(1) * SPESD 
SPEEDY = .8 * SPEEDY + NORM(2) * SPESD 
X = X + SPEEDX * DELTAT 
Y = Y + SPEEDY * DELTAT 
RANGE = SQRT(X**2 + Y**2) 
BEARI = ATAN2(Y,X) 
* GENERATION OF 2 NORMAL RANDOM BEARING AND RANGE 
CALL LNORPC(DSEED,NORM,2) 
MEASDV(1) = PERC * RANGE 
MEASDV(2) = .01 
IF(NOICE.EQ.1) THEN 
RANGEN = RANGE + NORM(1) * MEASDV(1) 
BEARIN = BEARI + NORM(2) * MEASDV(2) 
ELSE 
RANGEN = RANGE + ((NORM(1) * (1. + H * NORM(1) ** 2)) 
/ (QQ. + 6*H + 15* H**2)**.5)) * MEASDV(1) 
BEARIN = BEARI + ((NORM(2) * (1. + H * NORM(2) ** 2)) 
C / (1. + 6*H + 15* H**2)**.5)) * MEASDV(2) 
ENDIF 
XN = RANGEN * COS(BEARIN) 
YN = RANGEN * SIN(BEARIN) 
WRITE(10,115) TIME,X,Y, RANGE, 180.*BEARI/PI,XN,YN,RANGEN, 
C 180.*BEARIN/PI 
TIME = TIME + DELTAT 
390 CONTINUE 
GOTO 900 


* 


PEE SELL ELE TEE EE GENERATION OF A MODEL 4 TARGET SEEKHKEKRKSSKAEKERTKKEKERESE 
2 
400 WRITE(6,*) INPUT TARGET INITIAL X-Y POSITION (E.G., 0. 0.) 
READ(6,*) X,Y 
WRITE(6,*) INPUT MEAN TIME (HRS) BETWEEN COURSE CHANGES’ 
READ(6,*) ALPHA 
ALPHA = 1. / ALPHA 
WRITE(6,*) INPUT TARGET SPEED (KTS)’ 
READ(6,*) SPEED 
WRITE(6,*) INPUT RANGE STANDARD DEVIATION AS A PERCENTAGE OF” 
WRITE(6,*)) | ACTUAL RANGE (E.G. 5% IS 5) 
READ(6,*) PERC 
PERC = PERC / 100. 
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OPEN(UNIT =08, FILE=’D:\APL2\MOTIO4.RES’ , STATUS =’NEW’) 
* MAIN LOOP FOR MODEL 4 
DO 490 I=1,SIZE 
IF(ALTTIME .LE. TIME) THEN 
CALL LRNDPC(DSEED,U,1) 
COURSE = 2 * PI * U(1) 
CALL EXPOPC(DSEED,ALTT,ALPHA) 
ALTTIME = ALTTIME + ALTT 
SPEEDX = SPEED * COS(COURSE) 
SPEEDY = SPEED * SIN(COURSE) 
ENDIF 
X = X + SPEEDX * DELTAT 
Y = Y + SPEEDY * DELTAT 
TIME = TIME + DELTAT 
RANGE = SQRT(X**2 + Y**2) 
BEARI = ATAN2(Y,X) 
* GENERATION OF 2 NORMAL RANDOM BEARING AND RANGE 
CALL LNORPC(DSEED,NORM,2) 
MEASDV(1) = PERC * RANGE 
MEASDV(2) = .01 
IF(NOICE.EQ.1) THEN 
RANGEN = RANGE + NORM(1) * MEASDV(1) 
BEARIN = BEARI + NORM(2) * MEASDV(2) 
ELSE 
RANGEN = RANGE + ((NORM(1) * (1. + H * NORM(1) ** 2)) 
/ (1. + 6*H + 15* H**2)**.5)) * MEASDV(1) 
BEARIN = BEARI + ((NORM(2) * (1. + H * NORM(2) ** 2)) 
C / (1. + 6*H + 15* H**2)**.5)) * MEASDV(2) 
ENDIF 
XN = RANGEN * COS(BEARIN) 
YN = RANGEN * SIN(BEARIN) 
WRITE(08,115) TIME,X,Y,RANGE, 180.*BEARI/PI,XN, YN,RANGEN, 


Cc 180.*BEARIN/PI 
490 CONTINUE 
900 STOP 

END 


® 


* THIS SUBROUTINE WILL GENERATE A VECTOR OF NORMAL RANDOM 


* VARIABLES ACCORDING TO THE SINE-COSINE METHOD 
x 


SUBROUTINE LNORPC(DSEED,A,N) 
INTEGER N,],IND 

REAL A(N),U(2),S,W,U1,XSTAR 
DOUBLE PRECISION DSEED,PI 
DATA PI/3.14159265358979D0/ 

IND = 1 

DO 100 I=1,N 
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IND = -IND 
IF (IND.GE.0) GOTO 20 
10 CALL LRNDPC(DSEED,U,2) 
S =SQRT(-2*ALOG(U(1))) 
W = 2*PI*U(2) 
XSTAR = S*COS(W) 
A(I) = S*SIN(W) 
GOTO 100 
20 A(I) = XSTAR 
100 CONTINUE 
RETURN 
END 


SUBROUTINE LRNDPC (DSEED,U,N) 
INTEGER N, I 
REAL U(N) 
DOUBLE PRECISION D31M1, DSEED, D31 
C D31M1=2**31 - 1 
C D311 =2**31 
DATA D31M1/2147483647.D0/ 
DATA D31_ /2147483648.D0/ 
DO 5 I=1,N 
C DSEED = DMOQD(950706376.D0*DSEED,D31M1) 
DSEED = DMOD(16807.D0*DSEED,D31M1) 
5 Ud) = DSEED / D31 
RETURN 
END 
* 
* THIS SUBROUTINE GENERATES 1 EXPONENTIAL(LAMBDA) RANDOM 
* NUMBER 
SUBROUTINE EXPOPC(DSEED,EXPO,LAMBDA) 
DOUBLE PRECISION DSEED 
REAL EXPO,LAMBDA,UNI(1) 
CALL LRNDPC(DSEED,UNI,1) 
A = 1.0 - UNI(1) 
EXPO = - ALOG(A) / LAMBDA 
RETURN 
END 
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